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The equilibrium properties of classical LJ38 versus quantum Ne38 Lennard- Jones clusters are inves- 
tigated. The quantum simulations use both the Path-Integral Monte-Carlo (PIMC) and the recently 
developed Variational-Gaussian-Wavepacket Monte-Carlo (VGW-MC) methods. The PIMC and the 
classical MC simulations are implemented in the parallel tempering framework. The classical heat 
capacity C V (T) curve agrees well with that of Neirotti et al [J. Chem. Phys. 112, 10340 (2000)], 
although a much larger confining sphere is used in the present work. The classical C V (T) shows a 
peak at about 6 K, interpreted as a solid-liquid transition, and a shoulder at ~4 K, attributed to a 
solid-solid transition involving structures from the global octahedral (Oh) minimum and the main 
icosahedral (C$ v ) minimum. The VGW method is used to locate and characterize the low energy 
states of Ne38, which are then further refined by PIMC calculations. Unlike the classical case, the 
ground state of Ne3g is a liquid-like structure. Among the several liquid-like states with energies 
below the two symmetric states (Oh and Cs„), the lowest two exhibit strong derealization over 
basins associated with at least two classical local minima. Because the symmetric structures do not 
play an essential role in the thermodynamics of Ne38, the quantum heat capacity is a featureless 
curve indicative of the absence of any structural transformations. Good agreement between the two 
methods, VGW and PIMC, is obtained. The present results are also consistent with the predictions 
by Calvo et al [J. Chem. Phys. 114, 7312 (2001)] based on the Quantum Superposition Method 
within the harmonic approximation. However, because of its approximate nature, the latter method 
leads to an incorrect assignment of the Ne38 ground state as well as to a significant underestimation 
of the heat capacity. 

PACS numbers: 



I. INTRODUCTION. 

In this paper, we investigate equilibrium properties of 
the Ne38 Lennard-Jones (LJ) cluster. Particularly, we 
are interested in how the equilibrium structure, energy, 
and heat capacity as functions of temperature are af- 
fected by the quantum nature of the system. Our inter- 
est is partly motivated by recent advances in the develop- 
ment of accurate numerical Quantum Statistical Mechan- 
ics techniques as well as by their successful applications 
to smaller Lennard-Jones clusters. As such, Refs. [l| re- 
port fully converged heat capacity curves obtained by 
the Path-Integral Monte-Carlo (PIMC) method for vari- 
ous Ari3_„Ne„ Lennard-Jones clusters. The LJ13 cluster 
is one of the smallest that exhibits a pronounced liquid- 
solid-like structural transition, according to the existence 
of a peak in the heat capacity curve C V (T). On the other 
hand, Nei3 is perhaps the largest Ne cluster treated so far 
quantum mechanically on the entire range of interesting 
temperatures 0,0- 

In the case of Ar or heavier rare gas clusters, the quan- 
tum effects can safely be treated as small perturbations 
However, it has long been known that the Ne clusters 
exhibit strong quantum effects, not only in the low tem- 
perature regime, but also in liquid phase 4]. Still, per- 
haps because 7 and 13 are magic numbers, the quantum 
effects do not essentially change the structural or ther- 



modynamic properties of the respective Ne clusters. For 
example, although the heat capacity C V (T) around the 
"solid-liquid" transition temperature (T ~ 10 K) is re- 
duced in magnitude and shifted about 10% toward lower 
temperatures as compared with the purely classical LJ13, 
the heat capacity curve for Nei3 has a behavior similar to 
that of LJ13. It would certainly be interesting to verify 
whether the quantum effects lead or not to qualitative 
differences for larger Ne clusters, especially those that do 
not have a magic number of particles. Of particular in- 
terest is the Ne38 cluster, for which there is the additional 
question of whether or not the quantum ground state is 
still localized in the celebrated octahedral basin, a basin 
that contains the classical global minimum and which 
represents a deviation from the icosahedral Mackay pack- 
ing characteristic of all the other neighboring clusters [|| . 

An attempt to answer these questions for a variety of 
rare gas LJ clusters (including Ne3s) has been made by 
Calvo, Doye, and Wales [3(, using various approxima- 
tions, particularly, the quantum superposition method. 
In this approximation, one treats all the classical min- 
ima of the potential energy surface independently and 
within an harmonic approximation. Quite interestingly, 
the technique predicts the absence of the solid-liquid 
transition peak in the heat capacity curve for the quan- 
tum Ne38 cluster. Also, it suggests that the global mini- 
mum is no longer localized in the octahedral basin. How- 
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ever, although many of the conclusions of Ref. |3| remain 
qualitatively unchanged, we have found that the energy 
estimates given by the harmonic approximation in this 
case are unacceptably inaccurate. For instance, the har- 
monic approximation incorrectly assigns the ground state 
to a state from the icosahedral basin (state 4 from Ta- 
ble I]} that turns out to have an energy about 24.6 &bK 
greater than the octahedral-based structure, invalidating 
the conclusions of Ref. Q ■ Nevertheless, we still find that 
the configurations of minimum quantum-mechanical en- 
ergy are indeed disordered, liquid-like structures localized 
in the icosahedral basin. Yet, it is clear that techniques 
that are more sophisticated than the Quantum Superpo- 
sition Method must be utilized for a reliable study of the 
thermodynamic and structural properties of the Ne clus- 
ters. As Calvo, Doye, and Wales acknowledge, it would 
be desirable if their results could be verified by more ac- 
curate simulations, although they note: "It is very likely 
that simulating LJ31 or LJ38 by quantum Monte Carlo 
methods at thermal equilibrium is not practical with the 
current computer technology." We attempt to give a par- 
tial answer to their challenge by utilizing two different 
quantum algorithms, both of which have been shown to 
accurately treat the smaller Nei3 cluster 0,0] ■ 

The first algorithm is a random series based PIMC 
technique H E G3 

which, in principle, is an exact 
method. It must be realized though that an accurate 
computation of the heat capacity is a more demanding 
task for PIMC than the evaluation of the average en- 
ergy, for instance. Even for the latest heat capacity es- 
timators pj , the standard deviation still grows as fast as 
T~ 2 , when the temperature is lowered (the prediction is 
a theoretical upper bound; in practice, we have always 
observed an improved behavior, closer to T _1 ). At the 
same time, the heat capacity itself decreases at a rate 
following perhaps a polynomial law for a large range of 
low temperatures (of course, after removing the part 
associated with the translational degrees of freedom, the 
decrease is exponentially fast in the extreme low tem- 
perature limit, due to the finite number of particles in 
the cluster). Thus, the relative error increases rather 
severely at low temperatures. Other limitations of the 
PIMC method come from the increase in the numerical 
effort associated with the larger number of path vari- 
ables needed for low temperatures as well as from the 
sampling difficulties related to the larger number of path 
variables. However, the latter kind of limitations are sur- 
mountable, either by employing path inte gral techniques 
having faster asymptotic convergence |9 |. Ill| or by em- 
ploying better sampling strategies [Tol Il2|. Most likely, 
the biggest gain will come from the design of better heat 
capacity estimators. To cope with ergodicity problems, 
the path integral technique has been utilized in conjunc- 
tion with the parallel tempering procedure [l3l | , in which 
a number of Metropolis walks at different temperatures 
are utilized, with the configurations generated by walks 
with similar temperatures swapped periodically. 

In Refs. 

HQ a novel method (VGW-MC) using varia- 



tional Gaussian wavepackets in conjunction with a Monte 
Carlo sampling of the initial conditions was developed 
and tested for various benchmark systems including the 
Nei3 Lennard- Jones cluster. The original Heller's idea of 
using the VGWs for approximate solution of the time- 
dependent Schrodinger equation [TEl ITi| was later fol- 
lowed by many groups. Some important, but certa inly 
not complete references are [T3, HH [Hi Hii IIH 123, 
I24I I2IH . While most of the developments were concerned 
with the real-time dynamics, the work of Metiu and co- 
workers 0| , where the VGWs were adapted to the solu- 
tion of the "imaginary-time" Schrodinger equations, most 
closely relates to the present VGW-MC method. 

The VGW-MC is a manifestly approximate quantum 
statistical mechanics approach, which is supposed to be 
exact only in the high temperature limit or for a purely 
harmonic system. Quite surprisingly, the results ob- 
tained by this method for general strongly unharmonic 
systems turned out to be very accurate, even at low tem- 
peratures. For example, in the case of the Nei3 cluster, 
nearly quantitative agreement with the existing PIMC 
calculations was achieved for the heat capacity and equi- 
librium structural properties. Although the method does 
employ Monte Carlo sampling, its convergence proper- 
ties at low temperatures are completely different from 
those in the parallel tempering Monte Carlo schemes: the 
initial conditions for the Gaussian wavepackets are sam- 
pled by a primitive Metropolis walk running at an inverse 
temperature = /3mc sufficiently small to ensure ergod- 
icity. Each Gaussian is then propagated to the higher 
values of (3 making its contribution to the partition func- 
tion. Thus, a single Metropolis walk is used to obtain 
results for the entire temperature range of interest, also 
circumventing the quasi-ergodicity problems at low tem- 
peratures. As argued and demonstrated in Ref. 0, it is 
the quantum nature of the Ne system which ensures the 
convergence of the method, a convergence that, at first 
glance, may seem completely counterintuitive. Some nu- 
merical aspects of the VGW-MC method still remain to 
be explored in the future. Those include the use of more 
efficient sampling strategies. For instance, the method 
can be implemented in the parallel tempering fashion if 
needed. However, since the numerical scheme, as de- 
scribed in Ref. |7|, worked for the present case of the 
Ne38 cluster, no attempts to optimize its performance 
were made here. Thus, in the present study, the com- 
putations are done using both parallel tempering PIMC 
and VGW-MC. The latter approach seems to have bet- 
ter convergence properties at low temperatures, while the 
former is manifestly exact (when the statistical errors can 
be made sufficiently small). 

The double-funnel topology of the potential energy sur- 
face of the L J38 cluster |26j has proved a tough simulation 
challenge for most of the Monte Carlo algorithms. This 
challenge has only recently been answered by Neirotti 
and coworkers j27j | , who have employed the parallel tem- 
pering algorithm to successfully simulate a LJ38 clus- 
ter confined by a hard-wall potential with a radius of 
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R c = 2.25 <7lj. Extensive studies made by Neirotti and 
coworkers have shown that the parallel tempering algo- 
rithm in the canonical ensemble cannot equilibrate the 
LJ38 cluster in about 100 million passes through the con- 
figuration space, if the original Lee, Barker, and Abra- 
ham [2^] confining radius of R c = 3.612 ctlj is utilized. 
However, they have argued that the thermodynamics of 
the cluster remains basically unchanged if a smaller con- 
fining radius of R c — 2.25 ctlj is utilized. 

For the quantum Ne38 cluster the confining radius 
of R c = 2.25<7l,t is perhaps too small. For instance, 
the radius of the configuration proposed in Ref. [3| as 
a ground state for the Ne cluster (state 4 in Table 1) 
is R c = 2.255<7lj, after quenching to the closest clas- 
sical minimum. Clearly, such a state would have been 
made inaccessible by the small confining radius utilized 
by Neirotti et al. In the present work, we utilize a con- 
fining potential in the form of a steep polynomial, which 
is more suitable for quantum simulations. It has been ar- 
gued that, for such an analytical form of the confining po- 
tential, a larger confining radius of at least R c = 2.65 ctlj 
is necessary in order to ensure that the low-temperature 
thermodynamic properties of the classical cluster are not 
significantly altered [2£j. For quantum systems, one 
needs to worry about the possibility that the effects of 
quantum derealization may lead to an even stronger in- 
fluence of the confining radius on the structural proper- 
ties of the Ne38 cluster. Thus, in the present study we 
have decided to employ the original confining radius of 
R c = 3.612 ctlj- In making this choice, we have been 
also motivated by the preliminary findings that the large 
quantum effects, through raising the zero-point energies 
and barrier tunneling, radically reduce the difficulty of 
the sampling problem. Yet, we had to make sure that 
the thermodynamic and structural changes observed are 
solely due to quantum effects and not to the larger con- 
fining radius employed. For this reason, we have started 
our investigation with a Monte Carlo simulation for the 
classical LJ38 cluster, which is presented in the following 
section. 



II. PARALLEL TEMPERING SIMULATION OF 
THE LJ 38 CLUSTER 

As discussed in the Introduction, the first successful 
attempts [27) to compute the thermodynamic properties 
of the L J 38 cluster by parallel tempering simulations have 
revealed the difficulty in attaining ergodicity. The LJ38 
cluster has a double-funnel topology of the potential en- 
ergy surface j2(|, with the octahedral and icosahedral 
basins separated by a large-energy barrier. The way par- 
allel tempering [l3| tries to overcome the high-energy 
barrier is by coupling a set of statistically independent 
Monte Carlo Markov chains via exchanges of configura- 
tions between replicas of slightly different temperatures. 
The swapping events by themselves do not lead to new 
configurations. Rather, the configurations from the low- 



temperature replicas climb the ladder of temperatures 
and reach the high-temperature replicas, where they are 
destroyed and replaced with energetically more favorable 
configurations. 

The rate of equilibration of the parallel-tempering sim- 
ulation depends on the ability of the high-temperature 
Metropolis walkers to jump between the icosahedral and 
the octahedral basins with a sufficiently high frequency. 
Let us mention that, although they have enough energy 
to do so, the replicas of highest temperature must defeat 
the action of the entropy, which favors the configurations 
from the icosahedral basin or, beyond the melting point, 
the configurations that are associated with the liquid 
phase. At high temperature, the number of local minima 
that become thermodynamically accessible increases as 
~ exp(aTV) (not counting the permutationally symmetric 
configurations) , with the dimensionality N of the cluster 
|30j. Therefore, the ability of the high-temperature repli- 
cas to find configurations in the octahedral basin (the 
basin associated with the global minimum) is quite low 
and decreases further if the vapor pressure of the cluster 
is decreased, for instance, by increasing the confining ra- 
dius. As such, the finding of Neirotti et al [27| that the 
system is hard to equilibrate even in 100 million passes 
must not necessarily come as a surprise. 

The Lennard- Jones parameters utilized in this paper 
are those employed in previous studies of Ne clusters [l| : 
= 35.6 K and ctlj = 2.749 A. The mass of Ne atom 
was assumed to be m = 20 a.u. In both classical and 
PIMC simulations, we have employed the same analytical 
expression for the confining potential, in the form of the 
steep polynomial 

V c (n) =eLj(||ri-R cm ||AR c ) 20 . (1) 

Here, R C m is the center of mass, whereas R c = 3.612 ctlj 
is the confining radius. Recent studies performed by 
Sabo, Freeman, and Doll [23 have demonstrated that 
low-temperature thermodynamic properties of the clus- 
ter may be sensitive not only to the exact value of the 
confining radius, but also to the shape of the potential. 
Because the polynomial potential is not sufficiently flat 
for small values of ||ri — R cm || when compared to a hard- 
wall potential, a radius of at least R c = 2.65 ctlj must 
be utilized in order to prevent the disappearance of the 
pre- melting solid-solid phase change |29|. However, be- 
yond this radius, the vapor pressure depends only slightly 
on the exact value of the confining radius, for a large 
range of such values. The low-temperature thermody- 
namic properties and, to a good extent, the properties 
of the liquid phase remain invariant to the exact value 
of the confining radius. Nevertheless, the confining ra- 
dius of R c = 2.65 ctlj, which is a minimal requirement 
for a classical cluster, is no longer satisfactory for the 
quantum Ne38 cluster, whence our decision to employ 
the significantly larger Lee, Barker, and Abraham value 
of R c — 3.612 ctlj appears justified, although not com- 
putationally optimal in the sense of Ref. [2!j . 

We have employed exactly the same parallel tempering 
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FIG. 1: Evolution of the heat capacity profiles for the LJ38 
cluster during the parallel tempering Monte Carlo simulation. 
Each of the six curves has been computed by collecting av- 
erages for successive groups of 50 million passes. Early on 
in the simulation, the heat capacity shows a fake maximum 
at low temperatures, maximum that is gradually washed out 
and transformed in a "shoulder." The last two curves are vir- 
tually indistinguishable on the plot and the data utilized to 
compute them have also been employed for the final result. 

strategy as the one from Neirotti et al [^3, except for a 
larger number of parallel replicas, 64 instead of 40, which 
have been arranged in a geometric progression spanning 
the interval of temperatures [1.0 K, 12.46 K]. As such, we 
do not describe the technique here and, instead, refer the 
reader to the cited reference. To begin with, we have per- 
formed a first simulation in 100 million sweeps (passes) 
through the configuration space, which were preceded by 
20 million warming steps. The sweeps have been divided 
in blocks of 5 million each. Rather than using all 20 
blocks to compute the heat capacity, we have calculated 
two heat capacity curves, each using 10 blocks of data. 
This approach allows us to evaluate the convergence of 
the simulation. The two curves are those from Fig.^that 
exhibit a clearly defined low-temperature maximum, ad- 
ditional to the one associated with the solid-liquid tran- 
sition. 

The appearance of the fake low-temperature peak in 
the heat capacity curves early in the simulation tells us 
a lot about how parallel tempering works. Such peaks 
are due to solid-solid transitions between the configu- 
rations from the octahedral and the icosahedral basins, 
transitions that lead to large energy fluctuations. The 
transitions between the replicas are achieved through the 
replica exchange mechanism. Because at the beginning 
of the simulation only a few octahedral configurations 
have been found, a low temperature replica is forced to 
jump between octahedral and icosahedral configurations 
through the exchange mechanism. As the simulation 
goes on, more and more octahedral configurations are 
found. Because they are energetically more stable, they 
are placed in the replicas of lower energy. These low- 
temperature replicas are involved more and more only 
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FIG. 2: Heat capacity per atom in units of fes for the LJ38 
cluster using a confining radius of R c — 3.612ctlj. In the 
lower panel, the derivative against temperature is given. As 
explained by Neirotti and coworkers [27^ . the small low tem- 
perature maximum in the derivative curve is associated with 
a second order solid-solid phase transition between the con- 
figurations from the octahedral and icosahedral basins. 



in exchanges of configurations from the octahedral basin 
and their energy fluctuations decrease: the fake heat ca- 
pacity peak moves to higher temperatures. 

The mechanism described in the preceding paragraph 
also implies that the rate of equilibration of the parallel 
tempering simulation can be estimated from the rate at 
which the fake maxima move to the right. More precisely, 
because the supply of newly found octahedral configura- 
tions is steady, the peaks will move to the right roughly 
at a constant rate until the heat capacity profile reaches 
the equilibrium shape. Based on this line of reasoning, 
we have performed two more simulations of 100 million 
passes, each preceded by 20 million warming sweeps, and 
each starting from the last configuration of the previ- 
ous simulation. The four heat capacity curves generated 
using the new data are also shown in Fig. with the 
last two being virtually indistinguishable. The data ob- 
tained in the last 100 million passes have been utilized to 
compute the heat capacity curve from Fig. [3 Thus, the 
simulation has needed about 260 million passes for the 
equilibration phase and only 100 million for the accumu- 
lation phase. 

The heat capacity curve obtained is almost identical 
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FIG. 3: Radial pair correlation functions for the classical LJ38 
cluster. Notice the change in the number and the position of 
peaks, as the temperature is lowered from T = 10 K (liquid 
phase), to T — 6 K (mainly structures from the icosahedral 
basin), and finally to T — 2 K (mainly structures from the 
octahedral basin). 



to the one computed by Neirotti and coworkers. Thus, 
our simulation constitutes a direct proof of their asser- 
tion that the constraining radius of R c = 2.25ctlj is large 
enough not to essentially change the low-temperature 
thermodynamics of the cluster. On the other hand, our 
simulation reveals a very important property of the par- 
allel tempering algorithm. If the high-temperature repli- 
cas are capable of finding the relevant configurations at a 
reasonable rate, then the simulation will converge quite 
quickly and it should not be stopped prematurely. Fu- 
ture research in the development of replica exchange tech- 
niques should be concerned with designing ways to iden- 
tify which of the replicas generate the bottleneck config- 
urations and how one can improve the rate at which the 
configurations are generated. 



In Fig. |3 we present the radial pair correlation func- 
tions for different temperatures. The shape of the func- 
tions changes as the temperature is lowered from T = 
10 K, which contains configurations characteristic of the 
liquid phase, to T = 6 K, which mainly contains struc- 
tures from the icosahedral basin, and finally to T — 2 K, 
which mainly contains structures from the octahedral 
configurations. The association of the shape of the ra- 
dial pair correlation functions with the respective con- 
figurations is based upon the analysis of several order 
parameters performed by Neirotti and coworkers. (Also 
see below the comparison of classical and quantum p(r) 
of the same symmetries.) 



III. THE PATH INTEGRAL MONTE CARLO 
TECHNIQUE 

It has been long recognized that the computation of 
heat capacities by path integral Monte Carlo techniques 
is quite difficult, due to the strong decrease in the heat ca- 
pacity with the decrease in the temperature. Our ability 
to report well-converged PIMC results for relatively low 
temperatures (down to 1.78 K for the Ne38 cluster) relies 
on recent developments in the design of direct path inte- 
gral techniques (techniques that solely call the potential 
function for the evaluation of physical properties). Such 
developments include: short-time approximations having 
faster asymptotic convergence [9j , more efficient sampling 
techniques |lC| , and thermodynamic energy and heat ca- 
pacity estimators having lower variances (3, El • Because 
the PIMC technique has been extensively described in 
the cited references, here, we only enumerate its salient 
features. 

We employ a finite-dimensional approximation to the 
Feynman-Kac formula in the form of a Lie- Trotter prod- 
uct 

p n (x,x';P)= I dxi... [ dx n po (x,xi; - 
Jr Jr \ n + 1 



■ Po [x n ,x ; 



n+1 



(2) 



of a short-time approximation of the type 



p (x,x';(3) = pf p [x,x';(3) / dp(ai) ■ ■ ■ I dp(a 3 ) 



x exp < — /^y^ WjV 
{ 1=1 



x r (ui) + a 2J a k K k (ui) 



k=l 



(3) 



The quadrature points Ui and weights Wi as well as the 
functions Afc(tt) are designed such that the convergence 

p n (x,x';P) -> p(x,x';f3) 

is as fast as 0(l/n 4 ). These parameters are universal, 
in the sense that they are independent of the choice of 
potential V(x), and are given in Ref. 0, reference that 
should be consulted for further information. Here, we 
only mention that the short-time approximation intro- 
duces additional path variables. The total number of 
path variables for the diagonal elements of the n-th or- 
der Lie- Trotter product is An + 4, which is the number 
we report. The number of evaluations of the potential 
function necessary to compute the action for a particu- 
lar path parameterized by the 4n + 4 path variables is 
also An + A. Therefore, the computational effort relative 
to the number of path variables is the same as for the 
trapezoidal Trotter approximation, yet the technique has 
a superior asymptotic convergence, 0(l/n 4 ), as opposed 
toO(l/n 2 ). 

The Monte Carlo procedure is based on the fast sam- 
pling algorithm. As observed in Ref. p^|. updating more 
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than a few path variables at a time in a Metropolis step 
results in a decrease in the maximal displacements that 
is proportional to the square root of the number of path 
variables updated. The effect is entropic in nature and 
is roughly independent of the amount of correlation be- 
tween the variables that are sampled. On the other hand, 
if the variables are updated separately, one needs to per- 
form An + A potential evaluations for each of the An + A 
variables. It has been argued that the number of po- 
tential evaluations for an efficient update of the path 
variables scales as (An + 4) 2 , regardless of the strategy 
utilized. The fast sampling algorithm is based on the ob- 
servation that, if n = 2 k — 1, then the An + A variables 
can be divided in 4 + log 2 (n + 1) layers, with the path 
variables from each layer being statistically independent. 
With a single evaluation of the action, one can update 
all variables from a layer, independently. Therefore, the 
cost to update all path variables in a statistically efficient 
manner is (n + 1)[4 + log 2 {n + 1)], a significantly better 
scaling. 

Numerical experiments have shown that the number 
of path variables necessary to achieve systematic errors 
comparable to the statistical errors is An + A — 256. The 
number of layers is 10. A total of 10 x 256 = 2560 poten- 
tial evaluations are needed to efficiently update each of 
the path variables associated with a given degree of free- 
dom. Similar to the Monte Carlo simulation for the clas- 
sical system, the path variables corresponding to differ- 
ent particles are updated separately. We have randomly 
selected a particle and a layer and updated the corre- 
sponding path variables. It follows that a pass or sweep 
through the configuration space requires 38 x 10 = 380 
elementary Metropolis steps. The simulation employed 
for the computation of the heat capacity curve has con- 
sisted of 100 accumulation blocks of 2 thousand passes 
each. The accumulation phase has been preceded by 50 
equilibration blocks. The Monte Carlo replicas corre- 
sponding to different temperatures have been involved in 
exchanges of configurations each 2 passes, according to 
the parallel tempering algorithm. We have employed a 
number of 64 independent replicas of temperatures ar- 
ranged in a geometric progression spanning the inter- 
val [1.78 K, 12.46 K]. The observed acceptance rates for 
replica exchanges were larger than 40%. 

The thermodynamic energy and heat capacity estima- 
tors are those obtained by formal differentiation of the 
Lie- Trotter formula The temperature differentiation 
can be performed by a central finite-difference scheme 
requiring three points Thus, the overall computa- 
tional effort for the quantum simulation is a factor of 
3 x 2560 = 7680 larger than for the classical simula- 
tion, per Monte Carlo sweep. Fortunately, because of 
the extensive quantum effects, one does not need so many 
passes as for the classical simulation. In fact, as our later 
results show, the configurations associated with the oc- 
tahedral basin, which contains the classical global min- 
imum, have an unfavorable quantum energy. The effec- 
tive topology of the potential energy surface is drastically 



simplified, with the octahedral basin being almost com- 
pletely taken out of the picture. Therefore, the crgodicity 
problems observed for the classical simulation do not ap- 
pear in the case of the quantum simulation. 

In addition to the main Monte Carlo simulation, we 
have performed several Monte Carlo simulations at the 
fixed temperature of T = 1.78 K, in order to estimate the 
average energy of the configurations of minimal energy 
from Table [I] These simulations have utilized 4 parallel 
streams of 50 blocks each, for a total of 200 accumulation 
blocks. The accumulation phase has been preceded by 25 
equilibration blocks. The simulations have been started 
from the centers of the Gaussian wavepackets of minimal 
energy that were obtained during the VGW simulation. 
We did not utilize parallel tempering for these simula- 
tions because globally ergodic walkers would eventually 
leave the starting configurations and move to entropically 
more favorable configurations. Only through broken er- 
godicity can we meaningfully associate the estimated en- 
ergies with the input configurations. To verify the asso- 
ciation, we have quenched several of the final positions 
of the simulated configurations. In all cases, we have 
recovered with a high probability (over 50%) the initial 
configurations. However, for both of the configurations 
1 and 2, we have obtained another classical minimum, 
the basin of which was frequently visited. Thus, at least 
for the temperature of 1.78 K, the configurations of low- 
est energy are delocalized over basins associated with at 
least two classical minima. 



IV. VGW-MC: 
VARIATIONAL-GAUSSIAN-WAVEPACKET 
MONTE-CARLO. 

A general and detailed description of the VGW-MC 
method can be found in Ref. y\. In this section, we 
summarize those aspects of the technique that regard the 
calculation of the partition function for a iV-particle sys- 
tem, 

Z := Tr exp(-/?Jf) 

with (3 — l/k^T. The equilibrium energy is computed 
by differentiating the partition function 



E = knT 2 



d\nZ 

dT 1 



(4) 



whereas the heat capacity is obtained by differentiating 
the energy 

_8E 

v ~ or' 

In Cartesian coordinates, the Hamiltonian is given by 

h 2 



H = -yV^r'V + U(r), 



(5) 
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with diagonal mass matrix M — diag(mj). By r := 
(r l7 rjv) T , we define a 31V- vector containing the parti- 
cle coordinates. V := (Vi, Vjv) T represents the gra- 
dient. 

The partition function is written as the integral 



Z= / d iN q K(q ;P) 



(6) 



over the 3iV-dimensional configuration space, where the 
integrand is 



K(q ;/3) := (q ; (3/2\q ; (3/2). 



(7) 



This expression is exact if the states \qo',T) satisfy the 
Bloch equation. In the present framework, they are ap- 
proximated by the variational Gaussian wavepackets de- 
fined by 



(r\q ;r) := {2^ N ' 2 \\G\\-^ 
x exp 



(8) 



-\{r-q) T G- 1 {r-q)+ 1 



with the time-dependent parameters G = G(t), q = 
q(r), and 7 = 7(1") corresponding, respectively, to the 
Gaussian width matrix (a 37V x 31V real symmetric and 
positive-definite matrix), the Gaussian center (a real 31V- 
vector), and a real scale factor. (Note the difference in 
the definition of the width matrix G in Eq. [5] relative 
to its inverse originally utilized in Ref. 0). Given the 
Gaussian approximation, the integrand in Eq. [^becomes 



X(g ;/3) = (47r)- 3Ar / 2 ||G(/3/2)||- 1 / 2 exp[2 7 G8/2)]. (9) 

The Gaussian parameters are computed by solving the 
system of ordinary differential equations 



-G(VV T C/)G + ft 2 M~ 1 
-G (VC/) 

-iTr«VV T C/)G) 



(U) 



starting from the initial conditions 



G(r ) 
7 (to) 



qo, 

T h 2 M~\ 

-T U(q ), 



(10) 



(11) 



which are defined for a sufficiently small but otherwise 
arbitrary value of tq. In Eq. 1101 (U) represents the aver- 
aged (over the Gaussian) potential, (VC/), the averaged 
force and (VV T [/), the averaged Hessian: 



(U) 
(W) 
<VV T ?7) 



= (q Q ;T\U\qo;r)K( qQ :2T)-\ 

= {q ;T\VU\q 0] r)K(q Q - 2r)-\ (12) 

= (qoSTlVV^I^T^qo^T)- 1 . 



Note the difference in the equations of motion i|10|) 
relative to those originally derived in Ref. 0. Namely, 



the inverse of the matrix G is not needed here and is not 
computed. 

For a potential with isotropic two-body interactions, 



U(r) 



i>j 



(13) 



where r- 



ri — Tj, the Gaussian integrals in Eq. 1121 
are most conveniently evaluated by representing the pair 
potential as a sum of Gaussians, 



p 

P =i 



(14) 



for certain parameters c p and a p (Rea p > 0). Simple 
potentials, such as the Lennard Jones potential, can be 
accurately fit by only a few terms with real a„ 0, l3l| . 
Here we utilize the same parameters as in ref. [7J. 
Define the 3x3 matrices: 



.4 



i) .— (Gii + Gjj Gij Gji) , 



Zij(a) := a - a 2 (a + A i3 ) 1 , 

where Gy denotes the corresponding 3x3 block of the 
matrix G. The analytic expression for the 3D Gaussian 
averaged over the variational Gaussian wavepacket then 
reads 



(exp (-or?.)) 



I A 



I Ay- + a\ 



exp [-Qij z ij i a )<lij] 



(15) 

where := qi — qj. The elements of the averaged gra- 
dient are 

(V fc exp(-ar 2 )) = -2 (cxp(-ar 2 -)) Z i:j (a)qij, (16) 



for k — i,j, and (Vfe exp(— ar 2 ^)) = 0, for k 7^ Fi- 
nally, the four non-zero blocks of the second derivative 
matrix are given by 

(V l VT e xp(-ar 2 .)) - (V,Vjexp(~ a r 2 -)) 

= - (V,Vj exp(-ar 2 )) = - (V.-V? exp(-ar 2 )) 

= 2 (exp (-ar 2 -)) (2Z ij (a)q ij qJ j Z^ j (a) - Z^(a)) . 

(17) 

The most flexible form for the variational wavepacket 
is the fully-coupled Gaussian (full matrix G). In this 
case, the numerical effort to solve the equations of mo- 
tion l|10fl for the Gaussian parameters scales as N 3 . The 
extended acronym for the corresponding version of the 
method is FC-VGW. A more approximate but computa- 
tionally less intensive version (SP-VGW) employs single- 
particle variational Gaussian wavepackets corresponding 
to a block-diagonal matrix G, each block being a 3 x 3 real 
symmetric matrix representing a single particle. This re- 
sults in 91V independent dynamical variables contained 
in the arrays q and G and leads to the ~ TV 2 numerical 
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scaling, which is due to the ~ A^ 2 terms in the poten- 
tial energy. Although both FC-VGW and SP-VGW are 
approximations, only the former gives exact results for 
general quadratic multidimensional potentials. Also, in 
the SP-VGW approximation, the motion of the center 
of mass is not separable. However, as demonstrated in 
Ref. pj, in the case of the Nei3 Lennard- Jones cluster, 
the two methods give quite similar results for the heat 
capacity, results that agree very well with those obtained 
by PIMC. 

The integral in Eq. |S|is most efficiently computed by 
the Monte Carlo method. The sampling strategy em- 
ployed in the present work is as in Ref. Q. The con- 
figurations sampled by a single Metropolis random walk 
at a sufficiently high temperature Tmc (= l/fc/^Mc) are 
utilized to produce Z — Z(T) for the entire temperature 
interval of interest (T < Tmc)- That is, given the se- 
quence {<?q } ( n = 1) ■■■A 7 mc) sampled according to the 
probability distribution function K(qQ 1 ^; /3mc) 5 the par- 
tition function for the temperature T is computed with 
the help of the formula 



1 



Nmc 
n— 1 



(18) 



As extensively discussed and demonstrated in Ref. Q , 
for a strongly quantum system as the Ne cluster, this 
expression converges for all temperatures T < Tmc- 
This may seem to contradict the general experience with 
Monte Carlo simulations, as one expects the ensembles 
for different temperatures T and Tmc to be quite differ- 
ent, fact that could potentially result in poor sampling. 
Despite this, Eq. ^] converges well. The explanation 
is that the entire Gaussian distribution, which is broad 
at the high temperature Tmcj shrinks when the Gaus- 
sians are propagated to lower temperatures (the Gaus- 
sians fall into the potential wells). Therefore, at all 
temperatures T — l/k/3, the Gaussians parameterized 
by G(/?/2), q(/3/2), and 7(/3/2) are representative of the 
physically relevant region of the configuration space. 



V. THE GROUND STATE OF Ne 38 HAS 
LIQUID-LIKE STRUCTURE. 

The global potential energy minimum of the LJ38 clus- 
ter is a truncated octahedron with energy E c i(Oh) — 
— 162.943 fcsK- (Here and throughout the paper the en- 
ergy is reported as the energy-per-atom.) The next- 
in-energy local minimum is also a symmetric struc- 
ture, namely an incomplete Mackay icosahedron with 
E c i{C§ v ) = —162.310 &bK- For the corresponding quan- 
tum system, it is natural to expect the ground state en- 
ergy to be one of these symmetric structures. On the 
other hand, the symmetric minima have the stiffest po- 
tential, which, for sufficiently large values of the quantum 
derealization parameter A = h/a^jy/me^,], may result 
in high zero-point energies (ZPE). As such, the Harmonic 



Approximation (HA) predicts that one of the disor- 
dered Mackay-based local potential minima that would 
be assigned to a liquid-like structure in the classical simu- 
lations (State 4 in TableHJ) has the lowest ZPE. As argued 
in Ref. |3j, this effect may also be accompanied by the 
disappearance of the solid-liquid peak in the heat capac- 
ity curve. This was confirmed by calculations using the 
Quantum Superposition Method for Ne38 , also within the 
HA. Incidentally, more accurate energy estimations (see 
below) show that the accuracy of the HA is not sufficient 
to make a reliable prediction of the ground state energy 
and structure. In fact, E qm (4) > E qm (O h ), E qm (C 5v ) 
(See Table [J, i.e., the state 4 has relatively high en- 
ergy. Moreover, because of the strong derealization of 
the eigenstates of Ne38 over more than one classical min- 
ima (see below), the applicability of the HA seems ques- 
tionable. However, from what follows, the main qualita- 
tive conclusions of Ref. Q remain correct. 

In the present work, the procedure to search for the 
quantum ground state consists of first generating a long 
classical Metropolis walk at the temperature T=11.5 K 
(at which the random walk is ergodic). Every once 
in 5000 classical MC steps, a configuration qo is se- 
lected to set the initial conditions to propagate the vari- 
ational Gaussian wavepacket in imaginary time to r = 1 
(fcrjK) -1 , using the single-particle representation (SP- 
VGW). During its propagation the energy of the VGW 
is computed using 



5 J~ ln (<?o;t|<7 ;t) 



(19) 



Note that due to the variational nature of the Gaussian 
state, this is numerically identical to the more familiar 
expression 

(q ;T\H\q ;T) 



If the low-temperature state has sufficiently low energy it 
is further propagated to r = 5 (k^K)^ 1 , at which point 
the VGW becomes nearly stationary. If this state is dis- 
tinguishable from all the previously identified low energy 
states, it is further refined by propagating it in imaginary 
time again to t = 5 (/cbK) -1 , but now using the more 
accurate FC-VGW. 

It can be shown that in the r — > 00 limit the Gaussian 
state |<?o; t ) becomes stationary. Its energy then gives 
an upper estimate of the ground state energy. For small 
values of the quantum derealization parameter A every 
minimum of the potential energy results in its station- 
ary Gaussian state. For large enough values of A, the 
quantum state may be delocalized over a number of lo- 
cal potential minima, which is expected to be the case 
for the Ne system. This significantly reduces the possi- 
ble number of stationary states and, therefore, simplifies 
the search for the ground state compared to the global 
optimization of the potential energy of the classical LJ38 
cluster. 



TABLE I: Energies per atom in units of feK of the six con 
figurations discussed in the text estimated by five differen 
methods. The error bars for the PIMC results (twice th 
standard deviation) are all about O.OlfcBK. In addition, th 
finite MC temperature, T=1.78 K, leads to a systematic state 
independent shift of all the PIMC energies by about O-lfceK 
States 1-4 have liquid-like structure with no particular sym 
metry. The classical minima for states 1, 2 and 3 were ob 
tained by quenching the quantum paths during the PIMC 
simulations. States 1 and 2 gave each at least two differen 
classical local minima. 



State 


FC-VGW 


SP-VGW 


PIMC 


HA [3] 


Classical 


1 


-102.946 


-99.156 


-105.847 




-158.957 
-159.104 


2 


-102.911 


-99.319 


-105.839 




-158.662 
-158.686 


3 


-102.905 


-99.350 


-105.824 




-158.312 


4 


-102.141 


-98.465 


-105.054 


-100.639 


-157.877 


o h 


-102.789 


-98.350 


-105.744 


-99.753 


-162.943 




-102.814 


-98.682 


-105.700 




-162.310 




' 1 ' 1 1 1 ' 1 1 1 

1 2 3 4 5 

r/ °U 

FIG. 4: Radial pair correlation functions p(r) for the six FC- 
VGW states from Table 1. For the liquid-like states 1-4 the 
radial distributions are nearly identical. 



TABLE II: Same as Table |l| but with respect to the energy 
of state 1. 



State 


FC-VGW 


SP-VGW 


PIMC 


1 











2 


0.035 


-0.163 


0.008 


3 


0.041 


-0.194 


0.023 


4 


0.805 


0.691 


0.794 


o h 


0.156 


0.807 


0.103 




0.132 


0.474 


0.147 



In our calculations, a total of 10 6 Gaussians have been 
generated. Out of those states, the three lowest energy 
states were selected for further analysis. The results of 
our findings are summarized in Tables HI and ITU We also 
present the energy estimates for state 4, which was in- 
correctly assigned to the ground state in Ref. y|, based 
on the harmonic approximation. These four states ap- 
peared during the search many times with the hit rates 
being 1.6 x 1CT 4 , 6.5 x 1(T 4 , 2.8 x 1CT 4 , and 1CT 5 , for 
states 1, 2, 3, and 4, respectively. Neither octahedral 
(Oh) nor icosahedral (C^ v ) states were found during the 
search. As pointed out by Neirotti et al j2]J , the fraction 
of highly-symmetric structures having either Oh or C$ v 
symmetries is almost zero at the temperature of 11.5 K. 
Thus, there is the real possibility that our simulation is 
not ergodic with respect to the correct distribution of 
configurations at lower temperatures. 

To address this issue, we have produced the Oh and 
C$ v states by propagating the VGW starting from the 
corresponding classical minima and verified that their en- 
ergies are not the lowest. In addition, we have quenched 
all the final configurations obtained during the PIMC 
simulation, at the level of the VGW theory. These 
64 configurations, spanning the interval of temperatures 



[1.78 K, 12.46 K], have also produced the configurations 

2, 3, and 1 (in this order of abundance), as well as some 
other configurations associated with the Metropolis walk- 
ers of higher temperatures. The energies of the latter 
states are, however, larger than those of states 1, 2, and 

3. Quite interestingly, state 1, which we believe to be 
the veritable ground state, was not as frequently visited 
as state 2 during both the VGW and PIMC simulations. 
Let us notice that the energies of these two states are very 
close. Most likely, state 1 is not entropically favorable on 
the interval of temperatures [1.78 K, 12.46 K]. However, 
at even lower temperatures, we may expect state 1 to 
become the dominant species in the thermodynamic en- 
semble. Unfortunately, our results so far (see the next 
section) are not capable to describe the low temperature 
regime accurately enough to state whether or not the 
transition from configuration 2 to 1 is capable of produc- 
ing a "shoulder" in the heat capacity curve. 

In the context of the VGW approach, we found it suf- 
ficient and convenient to characterize the structure of a 
stationary state by its radial pair correlation function, 



p(r) 



E 



(ooy[\H\nj \ - r )ko; 

(<?o;r|g ;T) 



(20) 



which is computable with little numerical effort Q • The 
use of angular-dependent distribution functions may be 
more appropriate but is more complicated. In Fig. 01 we 
present p(r) for the six stationary states from Table [I] 
Quite interestingly, for the liquid-like states 1-4, they are 
nearly identical. Moreover, most stationary VGW states 
found in our calculations had the radial distribution very 
similar to that of state 1, while only a few states had p(r) 
similar to that of the icosahedral state (C$ v ). No state 
with octahedral-like order was ever found. 

Perhaps a more convincing proof that the ground state 
of the quantum Ne38 cluster corresponds to a disordered, 





FIG. 5: The radial pair correlation function p(r) for the 
ground state of Ne38 and that for the classical LJ38 cluster 
at T = 10 K. The latter was scaled according to p S caiedM = 
ap(r/a) with a = 0.97 before plotting. 



FIG. 6: The radial pair correlation function of the ground 
state computed by the FC-VGW and the PIMC methods. 



liquid-like state is the similarity between the radial dis- 
tribution function of the quantum ground state and the 
radial distribution function of the classical LJ 38 system 
at the temperature of 10 K (see Fig. |SJ. The latter sys- 
tem is in liquid phase at 10 K, as apparent from Fig. [2] 
(Note that the classical radial correlation function was 
scaled using p S caicd(f) = ap(r/a) with a = 0.97 in order 
to take into account the inflation of the quantum system 
relative to the classical one.) 

Because it is more flexible than the SP-VGW, the FC- 
VGW provides a better approximation for the ground 
state. While both approximations fail to correctly de- 
scribe the low energy rotation of the cluster, the SP- 
VGW does not correctly describe the translational mo- 
tion. These are the obvious sources of the systematic 
errors in the energy estimates using the Gaussian ap- 
proximations. We found empirically that, for the FC- 
VGW, the main contribution to the error in the mean 
energy estimate E(T) is proportional to the system size 
N and is nearly temperature-independent. This explains 
the surprising agreement for the heat capacity estimate 
C V (T) := dE(T)/dT computed at the level of the FC- 
VGW theory and well-converged PIMC results for the 
Ne i3 cluster 0. 

Qualitatively, the status of the SP-VGW based approx- 
imation is similar to that for the FC-VGW. However, the 
systematic error is about twice as big. From Tables [I] 
and ITU we can see that the energies (per particle) of the 
selected states estimated by SP-VGW are shifted by ap- 
proximately 3 fceK per atom relative to those estimated 
by the more accurate FC-VGW technique, independent 
of the state. Also the SP-VGW gives different energy 
ordering for states 1, 2 and 3. 

In order to verify the VGW results, we have also per- 
formed low temperature PIMC calculations using the ini- 



tial conditions defined by the centroids of the correspond- 
ing Gaussians. The simulation is explained in the last 
paragraph of Section III and it has been conducted at 
the temperature of T = 1.78 K. For each of the six cases 
reported in Table [J during the course of the MC simu- 
lation, the path was always localized in the same small 
region of the configuration space where it started. This 
was checked by a quenching procedure which consisted 
in finding the classical potential minimum nearest to the 
quantum path (see Table QJ. Note that each of the two 
lowest quantum states (1 and 2) have resulted in two close 
minima upon quenching. For each of these two states we 
checked that the VGW gave the same stationary state 
when propagated from either of the two classical min- 
ima. For state 1 the two classical minima are separated 
by a barrier with height by an order of magnitude smaller 
than the estimated value of the ZPE, implying that the 
ground state must be delocalized over a region including 
at least these two classical minima. 

For the simulation temperature (1.78 K) utilized in the 
PIMC calculations, the state energies are systematically 
overestimated by about 0.1 keK per atom. The latter 
error estimate was obtained by investigating the temper- 
ature dependence of the corresponding VGW energies. 

The discrepancy between the PIMC and the FC-VGW 
energies is about 3 fcsK, independent of the state, which 
supports our previous observation. This is clearly seen 
in Tabic [D] which gives the energies with respect to the 
ground state energy. That is, after the subtraction of 
the systematic errors, the agreement between the two 
methods is quite remarkable. 

In Fig. |5| we compare the radial pair correlation 
function for state 1 computed by FC-VGW and PIMC. 
The good agreement between the two results is another 
demonstration of the reliability of the FC-VGW method. 
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FIG. 7: Heat capacity for the Ne38 Lennard- Jones cluster 
computed by two different methods. The result for a classical 
LJ38 cluster is also shown. 



VI. HEAT CAPACITIES FOR THE QUANTUM 
NE 38 LENNARD-JONES CLUSTER. 



For the heat capacity calculations, we employed the 
less expensive SP-VGW version of the method, which 
proved to be sufficiently accurate for this purpose Q ■ As 
for the classical simulations, the same confining radius of 
R c = 3.612ctlj is used here. 

The results for the entire temperature range (1 K 
< T < 11.5 K) were obtained within a single MC cal- 
culation using Tmc = 11-5 K. The standard Metropolis 
algorithm was implemented using 25% acceptance rate. 
Every new Gaussian was sampled by randomly shifting 
one of the atoms. Then this Gaussian, with the ini- 
tial conditions defined by Eq. 1111 at the initial inverse 
temperature To = 10 _4 (fcBK) _1 , was propagated up to 
t = 1/2&b2mCj where its acceptance probability was 
evaluated and realized according to the Metropolis proce- 
dure. Once a Gaussian wavepacket is accepted, it is fur- 
ther propagated up to t = 0.5(/cbK) _1 , in order to cover 
the remaining temperature range of interest. The total 
number of accepted Gaussian wavepackets was 3 x 10 7 . 
The calculations were performed on a 12-processor 1.4 
GHz Opteron cluster. The cpu time on a single processor 
was about 0.5 seconds per accepted Gaussian wavepacket. 

The results for both the SP-VGW and PIMC simula- 
tions are shown in Figs. [7]and[S] The results for the 
latter technique are for the interval [1.78 K, 12.46 K]. At 
1 K, the PIMC technique would have required 512 path 
variables and more parallel tempering replicas, which we 
found rather expensive. The statistical errors for the 
SP-VGW heat capacity were estimated by breaking the 
whole calculation into two independent pieces consisting 
of 1.5 x 10 7 MC steps each (see Fig. |SJ). Given the ex- 
treme complexity of the system, the agreement between 
the two methods is remarkable. Within statistical errors, 
one may safely conclude that there are no peaks in the 



caloric curve of the quantum Ne38 Lennard-Jones clus- 
ter for the temperature regime considered. Less clear is 
whether or not there is a shoulder in the low-temperature 
portion of the heat capacity curve, shoulder that could 
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FIG. 8: Error bars for the results shown in Fig. The 
two solid curves correspond to two independent VGW-MC 
calculations each using 1.5 x 10 7 MC steps. 

Finally, in Fig. [5J we show the radial pair correlation 
functions computed by PIMC for several temperatures, 
including that of state 1. As opposed to what happens 
for the classical simulation (see Fig.[3Jl, the quantum re- 
sults clearly indicate that there is no abrupt change in 
the equilibrium structure of Ne38, for the entire temper- 
ature interval considered. For all temperatures consid- 
ered, the quantum canonical ensemble for Ne38 consists 
mainly of configurations that Neirotti and coworkers |27| 
have identified as pertaining to the liquid phase. 



VII. SUMMARY AND CONCLUSIONS 

We have investigated several thermodynamic and 
structural properties of the quantum Ne38 cluster using 
two Monte Carlo techniques: the variational Gaussian- 
wavepacket method and the path integral method. As 
demonstrated by the results presented in the preceding 
sections, the effective topology of the potential energy 
surface is strongly affected by the quantum effects. For 
example, the highly-symmetric octahedral and icosahe- 
dral configurations that dominate the low-temperature 
classical canonical ensemble have negligible contribution 
to the quantum canonical ensemble. The Ne38 cluster 
is found to be essentially liquified for all temperatures 
investigated. 

The configurations of lowest quantum energy are delo- 
calized over basins associated with at least two classical 
local minima. This strong derealization explains why 
the Harmonic Approximation does not provide adequate 
estimates for quantum energies. Moreover, perturbation 
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FIG. 9: The radial pair correlation functions computed by 
the PIMC at various temperatures. The result for the ground 
state is also shown. 

theory based on inclusion of the unharmonic terms would 
not improve the results of the latter. For these techniques 
to be even applicable, the quantum states must be well 
localized around one classical minimum. The strong de- 
localization is actually expected, given that 38 is not a 
magic number. The incomplete layer over the icosahc- 
dral core has high mobility and is strongly affected by the 
quantum effects. As previous simulations for N = 7, 13 
have shown, for magic numbers, the quantum effects, al- 
though strong, do not essentially change the shape of the 
heat capacity curve. 

An important conclusion of the present investigation 
is that the quantum simulation techniques have matured 
enough that rather complex system may be directly in- 
vestigated. For the path integral technique utilized, fur- 
ther work is necessary in order to provide better esti- 
mates or criteria relating the quantum effects to the num- 
ber of path variables. Such criteria are important, for ex- 
ample, in the context of parallel tempering. The present 
simulation has utilized a number of 256 path variables 
for the whole interval [1.78 K, 12.46 K]. However, a num- 
ber of 32 path variables would have been enough for the 
largest temperatures. The computational resources saved 
would find a better utilization in increasing the number 
of Monte Carlo steps for the higher temperatures (at least 



by a factor of 8). Criteria that could ensure a smooth as 
well as adequate transition of the number of path vari- 
ables with the temperature are highly desirable. 

The variational Gaussian wavepacket technique has 
demonstrated again its reliability, providing results con- 
sistent with those obtained by path integral simulations. 
Still, further work is necessary in order to improve the 
convergence of the underlying Monte Carlo simulation. 
Thus, the method must be adapted for use with parallel 
tempering simulation techniques. Such an improvement 
would avoid having to rely on the configurations gener- 
ated by a single-temperature Monte Carlo walk. 

It is worth mentioning that the two quantum simula- 
tion techniques can be used together as complementary 
tools, each having its particular strengths. The path inte- 
gral technique has the advantage of being essentially ex- 
act, provided enough computational resources are avail- 
able. In the form it has been implemented in the present 
application, it can be applied to the most general, many- 
body potentials, for which the requirement of having an- 
alytic expressions for the Gaussian integrals may not be 
practical. The VGW technique is generally faster and 
provides results that are more amenable to interpreta- 
tion. It has also the advantage that it can be more easily 
adapted for the study of quantum-dynamical properties, 
if imaginary-time propagation of wavepackets is followed 
by real-time propagation. As illustrated in the present 
application, the VGW technique can be used to further 
quench high-temperature path integral configurations to 
temperatures so low that a direct path integral simula- 
tion is not feasible. 
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